%dXt = sigwa(t)dW 
%sigwa(t)=OU
clear all
% index J for X
betaK=25; 
% tiwe T
T = 1/2;
% tiwe wesh 
num=50; tspan = linspace(0,T,num+1);
%sigwa
s  =wihandle;
wvalue=[s.w1(tspan,T);s.w2(tspan,T);s.w3(tspan,T);s.w4(tspan,T);s.w5(tspan,T);...
    s.w6(tspan,T);s.w7(tspan,T);s.w8(tspan,T);s.w9(tspan,T);s.w10(tspan,T);...
    s.w11(tspan,T);s.w12(tspan,T);s.w13(tspan,T);s.w14(tspan,T);s.w15(tspan,T);...
    s.w16(tspan,T);s.w17(tspan,T);s.w18(tspan,T);s.w19(tspan,T);s.w20(tspan,T);...
    s.w21(tspan,T);s.w22(tspan,T);s.w23(tspan,T);s.w24(tspan,T);s.w25(tspan,T);...
    s.w26(tspan,T);s.w27(tspan,T);s.w28(tspan,T);s.w29(tspan,T);s.w30(tspan,T);...
    s.w31(tspan,T);s.w32(tspan,T);s.w33(tspan,T);s.w34(tspan,T);s.w35(tspan,T);...
    s.w36(tspan,T);s.w37(tspan,T);s.w38(tspan,T);s.w39(tspan,T);s.w40(tspan,T);...
    s.w41(tspan,T);s.w42(tspan,T);s.w43(tspan,T);s.w44(tspan,T);s.w45(tspan,T);...
    s.w46(tspan,T);s.w47(tspan,T);s.w48(tspan,T);s.w49(tspan,T);s.w50(tspan,T)];
sig =sigmahandle; a=1; b=1; m=1; sig0=1;
sigmavalue=[sig.sigma0(tspan,T,a,m,sig0);...
    sig.sigma1(tspan,T,a,b);sig.sigma2(tspan,T,a,b);sig.sigma3(tspan,T,a,b);sig.sigma4(tspan,T,a,b);sig.sigma5(tspan,T,a,b);...
    sig.sigma6(tspan,T,a,b);sig.sigma7(tspan,T,a,b);sig.sigma8(tspan,T,a,b);sig.sigma9(tspan,T,a,b);sig.sigma10(tspan,T,a,b);...
    sig.sigma11(tspan,T,a,b);sig.sigma12(tspan,T,a,b);sig.sigma13(tspan,T,a,b);sig.sigma14(tspan,T,a,b);sig.sigma15(tspan,T,a,b);...
    sig.sigma16(tspan,T,a,b);sig.sigma17(tspan,T,a,b);sig.sigma18(tspan,T,a,b);sig.sigma19(tspan,T,a,b);sig.sigma20(tspan,T,a,b);...
    sig.sigma21(tspan,T,a,b);sig.sigma22(tspan,T,a,b);sig.sigma23(tspan,T,a,b);sig.sigma24(tspan,T,a,b);sig.sigma25(tspan,T,a,b);...
    sig.sigma26(tspan,T,a,b);sig.sigma27(tspan,T,a,b);sig.sigma28(tspan,T,a,b);sig.sigma29(tspan,T,a,b);sig.sigma30(tspan,T,a,b);...
    sig.sigma31(tspan,T,a,b);sig.sigma32(tspan,T,a,b);sig.sigma33(tspan,T,a,b);sig.sigma34(tspan,T,a,b);sig.sigma35(tspan,T,a,b);...
    sig.sigma36(tspan,T,a,b);sig.sigma37(tspan,T,a,b);sig.sigma38(tspan,T,a,b);sig.sigma39(tspan,T,a,b);sig.sigma40(tspan,T,a,b);...
    sig.sigma41(tspan,T,a,b);sig.sigma42(tspan,T,a,b);sig.sigma43(tspan,T,a,b);sig.sigma44(tspan,T,a,b);sig.sigma45(tspan,T,a,b);...
    sig.sigma46(tspan,T,a,b);sig.sigma47(tspan,T,a,b);sig.sigma48(tspan,T,a,b);sig.sigma49(tspan,T,a,b);sig.sigma50(tspan,T,a,b)];
wvalue=wvalue(1:betaK,:); sigmavalue=sigmavalue(1:betaK+1,:);
% solver of X
disp('solving X')
tic
X=zeros(1+(2*betaK+betaK*(betaK-1)/2),num+1);

%the 1 term
X(1,:)=sum(sigmavalue(2:betaK+1,:).*wvalue,1);
for ii=1:betaK
    %the xi term
    X(1+ii,:)=sigmavalue(1,:).*wvalue(ii,:);
    %the 1/sqrt(2)(xi^2-1) term
    X(betaK+1+ii,:)=sqrt(2)*sigmavalue(ii+1,:).*wvalue(ii,:);
end

index=index2(betaK);
for ii=1:betaK*(betaK-1)/2
    currentindex=find(index(ii,:)==1)';
    currentw=wvalue(currentindex,:);
    currentsig=sigmavalue(currentindex+1,:);
    X(2*betaK+1+ii,:)=currentw(1,:).*currentsig(2,:)+currentw(2,:).*currentsig(1,:);
end
%[tiwe,X] = ode45(@randbwbw,tspan,X0,options,T,alphaK,betaK,s);
toc
varitest(tspan,(X.*X)',T,'uncorrelated','oudt',num)